Suppression of thermal conductivity in graphene nanoribbons with rough edges 
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We analyze numerically the thermal conductivity of carbon nanoribbons with ideal and rough 
edges. We demonstrate that edge disorder can lead to a suppression of thermal conductivity by 
several orders of magnitude. This effect is associated with the edge-induced Anderson localization 
and suppression of the phonon transport, and it becomes more pronounced for longer nanoribbons 
and low temperatures. 



I. INTRODUCTION 

The study of remarkable properties of graphite struc- 
tures is one of the hot topics of nanoscience [1]. Graphene 
nanoribbons (GNRs) arc effectively low-dimensional 
structures similar to carbon nanotubes, but their main 
feature is the presence of edges. Due to the edges, 
graphene nanoribbons can demonstrate many novel prop- 
erties driven by their geometry, depending on their 
width and helicity. A majority of the current studies 
of graphene nanoribbons are devoted to the analysis of 
their electronic and magnetic properties modified by the 
presence of edges, including the existence of the localized 
edge modes [2, 3], which are an analog of surface states 
in the two-dimensional geometry. The edge can support 
localized vibrational states in both linear and nonlinear 
regimes [4, 5]. 

The effect of the edge disorder on the electronic trans- 
port of graphene nanoribbons has been discussed in sev- 
eral papers (see, e.g., Refs. [6-8]). It was found that al- 
ready very modest edge disorder is sufficient to induce the 
conduction energy gap in the otherwise metallic nanorib- 
bons and to lift any difference in the conductance be- 
tween nanoribbons of different edge geometry, suggest- 
ing that this type of disorder can be very important for 
altering other fundamental characteristics of GHRs. 

In addition to electronic properties, the thermal prop- 
erties of graphene are also of both fundamental and 
practical importance. Several experiments [9, 10] have 
demonstrated that graphene has a superior thermal con- 
ductivity, likely underlying the high thermal conductiv- 
ity known in carbon nanotubes [11]. This opens nu- 
merous possibilities for using graphene nanostructures in 
nanoscale thermal circuit management. 

Recent experiments demonstrated that thermal con- 
ductivity of silicon nanowircs can be dramatically re- 
duced by surface roughness [12, 13]. This results has 
been confirmed theoretically in the framework of a sim- 
plest phenomenological model of quasi-one-dimcnsional 
crystal that demonstrates the reduction of thermal con- 



ductivity due to roughness- induced disorder [14]. Molec- 
ular dynamics simulations [15] demonstrated that ther- 
mal conductivity of GRNs depends on the edge chirality 
and can be affected by defects. Therefore, we wonder 
if the edge disorder of GNRs can modify substantially 
their thermal conductivity, similar to the case of silicon 
nanowircs. 

In this Article, we study the thermal conductivity 
of isolated graphene nanoribbons with ideal and rough 
edges. By employing a direct modeling of heat trans- 
fer by means of the molecular-dynamics simulations, we 
demonstrate that the thermal conductivity grows mono- 
tonically with the GNR length as a power-law function. 
In contrast, rough edges of the nanoribbon can reduce 
the thermal conductivity by several orders of magnitude. 
This effect is enhanced for longer GNRs and for lower 
temperatures, and it corresponds to dramatic suppres- 
sion of phonon transport solely by the edge disorder. It 
means that nanoribbons with ideal edges can play a role 
of highly efficient conductors in nanocircuits, whereas the 
rough edges will transform them into efficient thermal re- 
sistors. 



II. MODEL 

We model a graphene nanoribbon as a planar strip 
of graphite, with the properties depending on the stripe 
width and chirality. The structure of the zigzag nanorib- 
bon can be presented as a longitudinal repetition of the 
elementary cell composed K atoms (the even number 
K > 4). Wc use atom numbering shown in Fig. 1(a). 
In this case, each carbon atoms has a two-component 
index a = (n, k), where n ~ 0, ±1, ±2, ... stands for the 
number of the elementary cells, and k = 1, 2, K stands 
for the number atoms in the cell. 

Each elementary cell of the zigzag nanoribbon has two 
edge atoms. In Fig. 1(a), we show these edge atoms 
as filled circles. We consider a hydrogen-terminated 
nanoribbon, where edge atoms correspond to the molecu- 
lar group CH. We consider such a group as a single effee- 
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FIG. 1: (a) Schematic view of a zigzag nanoribbon with 
rough edges and atom numbering. The edge atoms are shown 
as filled circles. Dotted lines separate the elementary cells of 
the nanoribbon. K is the number of atoms in the elementary 
cell, (b) Configurations of an ideal structure containing up to 
i-th nearest-neighbor interactions for i = 1, 5. 



tive particle at the location of the carbon atom. There- 
fore, in our model of graphene nanoribbon we take the 
mass of atoms inside the strip as Mq = I2m p , and for the 
edge atoms we consider a large mass Mi = 13m p (where 
m p = 1.6603 • 10~ 27 kg is the proton mass). 

To model two rough edges we randomly delete some 
atoms with second index k — 1 and k = K. Let < 
p < 1 be the probability of atom removal. As a result 
of the random atom removal from the edge layers, some 
atoms at the edges will have only one covalent bond C- 
C and should be deleted as well. After this operation, 
the edge become rough, as shown in Fig. 2. Here all 
edge atoms participate only in two valent bonds C-C. 
We characterize the degree of roughness by the parameter 
d = N a /Ni,, where N a is the number of atoms in an ideal 
nanoribbon, and iVj, is the number of atoms remaining in 
the edge-disordered nanoribbon after removing some of 
the edge atoms. Parameter d characterizes the density of 
the edge-disordered nanoribbon in comparison with the 
ideal case. When the probability p for removal of an 
edge atom is p = 0, we have d = 1, and d decays for 
larger values of the density p, so that for p = 1 (when all 
atoms with the second index k = 1, 2 and k = K — 1 , K 
are removed) , it takes the minimum value d = (K — 4) / K 
(for p = 1 we have again an ideal ribbon but for a smaller 
width, with K ~ 4 atoms in an elementary cell). For 
K = 12 and probability p = 0.5 the density is d = 0.87: 
this nanoribbon is shown in Fig. 2. 

To describe the dynamics of both ideal and disordered 
nanoribbons, we present the system Hamiltonian in the 
form, 
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where K — 4 < K„ < K is the number of atoms in the 
n-th elementary cell, M a is the mass of the hydrogen 



FIG. 2: Example of a zigzag nanoribbon with rough edges 
with N longitudinal segments. First left N+ segments are 
attached to the T = T+ thermostat and the last right N- 
segments are attached to the T = T_ thermostat. Number of 
atoms in the elementary cell K = 12, density of the nanorib- 
bon with rough edges d = 0.87 (probability p = 0.5). 



atom with the index a = (n, k) (for internal atoms we 
take M a = Mq, whereas for the edge atoms we take a 
larger mass, M a = Mi > Mo), u Q = (x a (t),y a (t), z a (tj) 
is the radius- vector of the carbon atom with the index a 
at the moment t. The term P a describes the interaction 
of the atom with the index a = (n, k) with its neigh- 
boring atoms. The potential depends on variations of 
bond length, bond angles, and dihedral angles between 
the planes formed by three neighboring carbon atoms, 
and it can be written in the form 

p = E f/ i + E c/ 2+E^+E c/ 4 + E c/ ^ (2) 

where f2j, with i = 1,2,3,4,5 stand for the sets of con- 
figurations including up to nearest-neighbor interactions. 
Owing to a large redundancy, the sets only need to con- 
tain configurations of the atoms shown in Fig. 1(b), in- 
cluding their rotated and mirrored versions. 

The potential U\ (u a , ) describes the deformation en- 
ergy due to a direct interaction between pairs of atoms 
with the indices a and /3, as shown in Fig. 1(b). The po- 
tential [^(iia, up, u 7 ) describes the deformation energy 
of the angle between the valent bonds u Q u^ and u,gu 7 . 
Potentials Ui(u a , up, u 7 , us), i = 3, 4, 5, describes the 
deformation energy associated with a change of the effec- 
tive angle between the planes u Q , u^, u 7 and u^, u 7 , u^. 

We use the potentials employed in the modeling of the 
dynamics of large polymer macromolecules [16, 17] for 
the valent bond coupling, 

U 1 {u 1 ,u 2 ) = ei{exp[-ao(/>-/3o)] - I} 2 , P = |u 2 -ui|, 

(3) 

where t\ = 4.9632 eV is the energy of the valent bond 
and po = 1.418 A is the equilibrium length of the bond; 
the potential of the valent angle 

J7 2 (ui,u 2 ,u 3 ) = e 2 (cos(^ - cos<y9 ) 2 , (4) 
cos<^ = (u 3 - u 2 ,ui - u 2 )/(|u 3 - u 2 | • |u 2 - Ui|), 

so that the equilibrium value of the angle is defined as 
cosv?o = cos(27r/3) = —1/2; the potential of the torsion 
angle 



J7 l (ui,U 2 ,U 3 ,U 4 ) = 6i(l - z,-cos< 



(5) 
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cos0 = (vi,v 2 )/(|vi| • |v 2 |), 
Vl = (u 2 - Ui) x (u 3 - u 2 ), 
V2 = (u 3 - u 2 ) x (u 3 - u 4 ), 

where the sign z,; = 1 for the indices i = 3, 4 (equilibrium 
value of the torsional angle <f>o = 0) and Zi = — 1 for the 
index i = 5 (0o = tt). 

The specific values of the parameters are ao — 
1.7889 A" 1 , e 2 = 1-3143 eV, and e 3 = 0.499 cV, and 
they are found from the frequency spectrum of small- 
amplitude oscillations of a sheet of graphite [18]. Ac- 
cording to the results of Ref. [19] the energy €4 is close to 
the energy e 3 , whereas £5 <C €4 (|e 5 /e4| < 1/20). There- 
fore, in what follows we use the values 64 = e 3 = 0.499 
eV and assume 65 = 0, the latter means that we omit the 
last term in the sum (2). 

III. METHODS 

In order to model the heat transport, we consider 
the nanoribbon of a finite length with two ends places 
in thermostats kept at different temperatures, as shown 
schematically in Fig. 2. In order to calculate numerically 
the coefficient of thermal conductivity, we should calcu- 
late the heat flux at any cross-section of the nanoribbon. 
Therefore, first we obtain the formula for calculating the 
longitudinal local heat flux. 

We define the 3A„-dimensional coordinate vector u„ = 
{xn,k, Un,k, z n,k}k=i which determines the atom coordi- 
nates of an elementary cell n, and then write the Hamil- 
tonian (1) in the form, 

H = y^fen = ^[-(M n Un,u„) + P„(u n _i,u„,u n+ i)], 

71 71 

(6) 

where the first term describes the kinetic energy of the 
atoms (M„ is diagonal mass matrix of the n-th elemen- 
tary cell), and the second term describes the interaction 
between the atoms in the cell and with the atoms of 
neighboring cells. 

Hamiltonian (6) generates the system of equations of 
motion, 

— M„U„ = F„ = Pi : „ + i + P2,n + P3,n-1; (7) 

where the function P 4i „ = Pj(u n _i, u„, u n+ i), Pi = 
<9P(ui, u 2 , u 3 )/<9ui, % = 1, 2, 3. 

Local heat flux through the n-th cross-section, j n , de- 
termines a local longitudinal energy density h n by means 
of a discrete continuity equation, h n = j n — jn+i- Using 
the energy density from Eq. (6) and the motion equa- 
tions (7), we obtain the general expression for the en- 
ergy flux through the n-th cross-section of the nanotube, 

jn = (Pl,n,U„_i) — (P 3i „_i,U„). 

For a direct modeling of the heat transfer along the 
nanoribbon, we consider a nanoribbon of a fixed length 
(N — l)h with fixed ends. We place the first N + = 40 



segments into the Langevin thermostat at T + = 31 OK, 
and the last iV_ = 40 segments, into the thermostat at 
T_=290K - see Fig. 2. As a result, for modeling of the 
thermal conductivity we need integrating numerically the 
following system of equations, 

M„u„ = -F„-rM„u„ + S+ for n = 2,...,N + , 
M n u n = -F„, for n = JV + + l,...,JV-JV_, (8) 
M„u„ = -F„-rM„u„ + H-, 

for n = N -N- +1,...,N, 

where T = l/t r is the damping coefficient (relaxation 
time t r = 0.1 ps), and 

is 12A"„-dimensional vector of normally distributed ran- 
dom forces normalized by conditions 

<^i(*i)C£(*a)) = 2M n , i k B T±8 nl 8 ij 5(t 1 - t 2 ). 

Details of the numerical procedure for modeling of ther- 
mal systems can be found elsewhere [20]). 

We select the initial conditions for system (8) corre- 
sponding to the ground state of the nanoribbon, and solve 
the equations of motion numerically tracing the transi- 
tion to the regime with a stationary heat flux. At the 
inner part of the nanotube (N + < n < N — iV_), we 
observe the formation of a temperature gradient corre- 
sponding to a constant flux. Distribution of the average 
values of temperature and heat flux along the nanotube 
can be found in the form, 

1 f* 

T n = hm / (M„u„(r),u„(r))dr, 

t^oo AK n k B t J 

h f* 

J n = lim - / j n (r)d,T, 

t->oo t J 

where fcs is the Boltzmann constant. For nanoribbons 
with rough edges we make the averaging not only in time 
but also on 240 independent realizations of the roughness. 

Distribution of the temperature and local heat flux 
along the rough-edged nanoribbon is shown in Figs. 
3(a,b) and 4(a,b). The heat flux in each cross-section 
of the inner part of the nanoribbon should remain con- 
stant, namely J n = J for N + < n < N — N_. The re- 
quirement of independence of the heat flux J„ on a local 
position n is a good criterion for the accuracy of numer- 
ical simulations, as well as it may be used to determine 
the integration time for calculating the mean values of J n 
and T n . As follows from the figures, the heat flux remains 
constant along the central inner part of the nanoribbon. 

A linear temperature gradient can be used to define 
the local coefficient of thermal conductivity, k(N — N + — 
N-) = {N-N--N+- 1)J/(T N++1 -T N - N _)S, where 
S = 2(dD y + 2rc)fc is the area of the nanoribbon cross- 
section (nanoribbon width D y = (3A/4 — l)po, Van der 
Waals carbon radius r c = 1.85A). Using this definition, 
we can calculate the asymptotic value of the coefficient 
k = limAr^oo k(N). 
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FIG. 3: (Color online) Distribution of (a) local heat flux 
J n and (b) local average temperature T n along ideal zigzag 
nanoribbon with K = 8. Length of the nanoribbon is 
L = (TV - l)h = 58.8 nm (TV = 240, h = 0.246 nm), and 
temperatures are T+ = 310 K and T_ = 290 K, the numbers 
of end segments interacting with the thermostats N± — 40 
(corresponding fragments are shown in grey). Heat conduc- 
tivity is k = 177 W/mK. 



IV. RESULTS 

Our analysis of linear eigenmodes of the nanoribbon 
with periodic boundary conditions in n reveals that in 
the case of edge disorder almost all vibrational modes 
are localized as functions of the longitudinal index n. 
This means that in our system we observe the manifes- 
tation of the Anderson localization due to the edge dis- 
order, earlier discussed only for the wave transmission in 
surface-disordered waveguides [21, 22]. 

To analyze oscillation eigenmodes, we define the dis- 
tribution function of the oscillatory energy along the 
nanoribbon as follows: 



E 

fe=r 



M nifc (|e„, M | 2 + |e„. fe . 2 | 2 + K M | 2 )/M , 



where n = 1,2, ...,N, M n ^ is mass of the atom with 
the index (n,k), and {e ra ,fe,i}f = i is a component of the 
corresponding eigenvector (see Ref. [4]). The energy 
distribution is normalized in accord with the condition: 
J2n=iPn = I- To describe the longitudinal energy local- 
ization, we introduce a new parameter, D ~ 1/ Xm=i P„, 
that characterizes the width of the energy localization 



FIG. 4: (Color online) Distribution of (a) local heat flux J n 
and (b) local average temperature T n along a zigzag nanorib- 
bon with rough edges (K — 10, density d — 0.87). Length 
of the nanoribbon is L = (TV - l)h = 58.8 nm (TV = 240, 
h — 0.246 nm), and temperatures are T+ = 310 K and 
TL = 290 K. Heat conductivity is k = 14 W/mK. 



along the nanoribbon. If a vibrational mode is localized 
only on one elementary cell, the corresponding width is 
D = 1. In the opposite limit, when the vibrational en- 
ergy is distributed equally on all elementary cells, we have 
D = TV, so that in a general case 1 < D < TV. 

Dependence of the width D on the frequency of the 
oscillatory eigenmodes is shown in Fig. 5. For an ideal 
nanoribbon, all modes are not localized: when the length 
TV = 300 we have the width 200 < D < 300. For nanorib- 
bons with rough edges, only the modes with the wave- 
length of the order of the nanoribbon length are not lo- 
calized. As a result, we expect that the edge disorder 
should lead to suppression of phonon transport and dra- 
matic reduction of the thermal conductivity. 

Our numerical results demonstrate that the thermal 
conductivity of graphene nanoribbon depends crucially 
on the degree of edge roughness. In spite of the fact that 
the nanoribbon has an ideal internal structure, its ther- 
mal conductivity is reduced dramatically, and it becomes 
much lower that the conductivity of an ideal nanoribbon 
of the same width. 

Distribution of the thermal flow J„ and local temper- 
ature T n along the ideal nanoribbon and the nanoribbon 
with rough edges (for the density d = 0.87) are presented 
in Figs. 3(a,b) and 4(a,b). In comparison with the ideal 
nanoribbon, the edge disorder leads to reduction of the 
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FIG. 5: (Color online) Dependence of the width parameter D 
of linear eigenmodes for the nanoribbon with periodic bound- 
ary conditions (rj + JV = n, length N = 300) on the frequency 
uj for an ideal nanoribbon (curve 1, width K = 10) and a 
nanoribbon with rough edges (curve 2, K = 10, p = 0.5). 
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FIG. 6: Dependence of the coefficient of thermal conductivity 
k of a finite nanoribbon with rough edges (N = 240, N± — 40, 
K = 12) on the density d. 



thermal flow in at least ten times, as well as it changes the 
temperature profile along the nanoribbon. In addition, 
in an ideal nanoribbon we observe thermal resistance at 
the edges placed into a thermostat, which disappears in 
the case of rough surfaces. As a result, for the length 
L = (N -N-- N + )h = 39.4 nm, (N = 240, N± = 40) 
the coefficient of thermal conductivity of the nanoribbon 
with rough edges is found as K = 14 W/mK that is in 
12.6 times lower than the thermal conductivity of an ideal 
nanoribbon, k = 177 W/mK. 

Dependence of the coefficient of thermal conductivity 
k on the degree of roughness characterized by the pa- 
rameter d is shown in Fig. 6 for K = 12, N = 240, 
and N± = 20. As follows from this figure, the ther- 
mal conductivity will be the lowest for the densities 
0.76 < d < 0.93 (corresponding to the probability of 
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FIG. 7: Dependence of the coefficient of thermal conductiv- 
ity n on the length of the central part of the nanoribbon L 
(dimension [L] =A). Curve 1 corresponds to an ideal nanorib- 
bon with K = 8, curve 2 - to edge-disordered nanoribbon with 
K = 10, p= 0.5. 



removing the edge atoms, 0.2 < p < 0.7). The max- 
imum is observed for d = 1 (probability p = 0) and 
d = (K — A)/K = 2/3 (p = 1) when we have ideal 
nanoribbons with K = 12 and K = 8 atoms in an el- 
ementary cell, respectively. The minimum is observed 
for the density d = 0.84 (probability p = 0.4). Below, 
we consider nanoribbons with rough edges created by re- 
moving edge atoms with the probability p = 0.5. The 
corresponding structure of this nanoribbon is shown in 
Fig. 2. 

Our numerical modeling described above demonstrates 
that for T = 300K the thermal conductivity of an ideal 
nanoribbon grows with its length L as a power-law func- 
tion, k ~ L a for L — > oo where a w 1/3. In contrast, 
the thermal conductivity of a nanoribbon with rough 
edges grows much slower, see Fig. 7. This difference 
grows with the length of the nanoribbon. For exam- 
ple, for L = 4.91 nm (N = 120, N ± = 40), a ratio /3 
between the coefficient of thermal conductivity of dis- 
ordered (K = 12, p = 0.5) nanoribbon, K\, and ideal 
(K = 10, p = 0) nanoribbon, Kq, of the same width is 
j3 — Ki/kq = 0.14, but for the length L = 314.4 nm this 
ratio becomes much smaller, j3 = 0.06. 

Efficiency of the thermal conductivity of edge- 
disordered nanoribbons also decreases with temperature, 
as well as with the ratio j3. For an ideal nanoribbon, 
the coefficient of thermal conductivity grows monotoni- 
cally for low temperatures (see Fig. 8, curve 1), so that 
for T — > we obtain k — > oo. This is related to the 
fact that the dynamics of nanoribbons approached the 
dynamics of one-dimensional linear system with infinite 
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FIG. 8: Dependence of the thermal conductivity coefficient 
k on the temperature T for the ideal finite nanoribbon (curve 
f, width K = 8) and for edge-disordered nanoribbon (curve 
2, width K = 10, probability p = 0.5). Length of the central 
part L = 320/i = 78.6nm (TV = 400, N± = 40). 



ductivity decays, since a diffusion transport is driven by 
nonlinear dynamics. As a result, the ratio (3 = K\/ko 
decays monotonically. For example, for the nanoribbon 
with the length L = 78.6 nm, at T = 500K this ratio is 
(3 = 0.10, at T = 300K, we have f3 = 0.077, at T = 100K 
we obtain [3 = 0.036, at T = 50K we find /3 = 0.021, 
and at T = 25K, we have j3 = 0.012 (i.e. the thermal 
conductivity is reduced by two orders!). 



V. CONCLUSIONS 

We have studied numerically thermal conductivity of 
carbon nanoribbons with ideal and rough edges. We 
have demonstrated that thermal conductivity of an ideal 
nanoribbon is a monotonic power-like function of its 
length. However, the thermal conductivity is modified 
dramatically when the structure of nanoribbon edges 
change. In particular, the thermal conductivity of a 
nanoribbon with an edge-induced disorder is reduced by 
several orders of magnitude, and this effect is more pro- 
nounced for longer ribbons and low temperatures. As 
a result, nanoribbons with ideal edges can play a role of 
highly efficient conductors, while nanoribbons with rough 
edges become efficient thermal resistors. 



thermal conductivity. In contrast, for the nanoribbon 
with rough edges we observe that for T > 100K its ther- 
mal conductivity depends only weakly on temperature, 
see Fig. 8, curve 2. This result is explained by the fact 
that in the edge-disordered nanoribbon all linear vibra- 
tional modes becomes localized due to the edge disorder, 
and the phonon transport is suppressed. For low tem- 
peratures, the system become linear and its thermal con- 
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